######### Figures in main text ########

### Figure 1 austerity_6m * dunemployment ###
attach(cabinet_filter_data)
set.seed(1234567)

int.mod.2<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_6m*d.unemployment,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.2)



## Instantaneous effect ##

pred_effects_d.dunemployment <- function(d.dunemployment){
  betas <- coef(int.mod.2)
  betas[1,2] <- 0.6268281
  betas<-as.numeric(betas[c(1), ])   
  varcov <- vcov(int.mod.2)
  set.seed(1234567)
  beta.tilde <- data.frame(mvrnorm(1000,betas,varcov))
  beta.tilde <- beta.tilde[,c(2,3,8)]
  
  colnames(beta.tilde) <- c(
    "alpha_1","beta_0","beta_1") 
  
  # create relevant comparisons 
  dunemployment <- seq(unname(summary(cabinet_filter_data$d.unemployment)[1]),
                       unname(summary(cabinet_filter_data$d.unemployment)[6]),length.out=250) # vary level
  # create hypothetical values based on these
  hyp_zt <- dunemployment
  hyp_zt_l1 <- dunemployment - d.dunemployment
  hyp_zt_p1 <- dunemployment + d.dunemployment
  # create list to loop over
  hyplist <- as.list(NULL)
  # create storage for results
  effects <- data.frame(matrix(ncol=7,nrow=length(dunemployment)))
  colnames(effects) <- c("dunemployment","inst.lo","inst.est","inst.hi","total.lo","total.est","total.hi")
  for(i in 1:length(dunemployment)){
    ### create temp storage
    hyplist[[i]] <- beta.tilde
    
    ### create hypotheticals
    # lagged level
    
    hyplist[[i]]$beta_1_hyp_zt_l1 <- hyplist[[i]]$beta_1*hyp_zt_l1[i]

    # current level
    hyplist[[i]]$beta_1_hyp_zt <- hyplist[[i]]$beta_1*hyp_zt[i]

    # future level
    hyplist[[i]]$beta_1_hyp_zt_p1 <- hyplist[[i]]$beta_1*hyp_zt_p1[i]

    # negative future level
    hyplist[[i]]$beta_1_hyp_zt_p1_neg <- (-1)*hyplist[[i]]$beta_1*hyp_zt_p1[i]
    
    ### calculate quantities of interest
    hyplist[[i]]$inst.effect <- rowSums(hyplist[[i]][,c("beta_0", "beta_1_hyp_zt")])
    
    hyplist[[i]]$total.effect <- (rowSums(hyplist[[i]][,c("beta_0", "beta_1_hyp_zt"
    )]))/(1-(hyplist[[i]]$alpha_1))
    effects$dunemployment[i] <- hyp_zt[i]
    effects$inst.lo[i] <- unname(quantile(hyplist[[i]]$inst.effect, .025))
    effects$inst.est[i] <- unname(summary(hyplist[[i]]$inst.effect)[4])
    effects$inst.hi[i] <- unname(quantile(hyplist[[i]]$inst.effect, .975))
    effects$total.lo[i] <- unname(quantile(hyplist[[i]]$total.effect, .025))
    effects$total.est[i] <- unname(summary(hyplist[[i]]$total.effect)[4])
    effects$total.hi[i] <- unname(quantile(hyplist[[i]]$total.effect, .975))
  }
  return(effects)
}

effects <- pred_effects_d.dunemployment(d.dunemployment = 0)

require(dplyr)
distribution.dta<-select(cabinet_filter_data,c(d.unemployment,protest_freq,imf))

ins.dunemployment.austerity <- ggplot(effects, aes(x=dunemployment,y=inst.est)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=inst.lo, ymax=inst.hi), alpha=.15) +
  geom_line(aes(x=dunemployment,y=inst.est)) +
  coord_cartesian(ylim=c(-5,2.5), xlim = c(-5,5)) + 
  geom_hline(yintercept=0, linetype=2) + 
  labs(title="Instantaneous effect of austerity_6m", x=expression(paste(Delta, "unemployment")),y="Vote Intention (cabinet)" 
       +theme(plot.title = element_text(size=24),axis.title= element_text(size=20))) +
  geom_rug (data=distribution.dta,aes(x = d.unemployment),
            inherit.aes = F, alpha = .5,
            position = "identity",size=1,sides="b")
ins.dunemployment.austerity



ins.dunemployment.austerity <- ggplot(effects, aes(x=dunemployment,y=inst.est)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=inst.lo, ymax=inst.hi), alpha=.15) +
  geom_line(aes(x=dunemployment,y=inst.est)) +
  coord_cartesian(ylim=c(-5,2.5), xlim = c(-5,5)) + 
  geom_hline(yintercept=0, linetype=2) + 
  labs(title="Instantaneous effect of austerity_6m", x=expression(paste(Delta, "unemployment")),y="Vote Intention (cabinet)") +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20)) +
  geom_histogram (data=distribution.dta,aes(x = d.unemployment,y = (-..count..)*0.07),
                  inherit.aes = F, alpha = .3, fill="grey",
                  position = "stack", bins=30,binwidth = 0.1)
ins.dunemployment.austerity


## Temporal illustration ## 

cabinet_filter_data[,"austerity_dunemployment"]<-cabinet_filter_data$austerity_6m*cabinet_filter_data$d.unemployment
attach(cabinet_filter_data)
set.seed(1234567)
int.mod.2<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_dunemployment,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.2)

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.25, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.5, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

Scen3 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = quantile(cabinet_filter_data$d.unemployment,probs =0.75, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(protest_freq,
                                      na.rm = TRUE),
                    austerity_dunemployment=0)

shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_dunemployment<- 0
shock_scen1$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.25, na.rm = T)

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_dunemployment<- 0
shock_scen2$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.5, na.rm = T)


shock_scen3 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen3$austerity_6m [6:11]<- 1 
shock_scen3$austerity_dunemployment<- 0
shock_scen3$austerity_dunemployment[6:11]<-1*quantile(cabinet_filter_data$d.unemployment,probs =0.75, na.rm = T)

set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
set.seed(1234567)
Sim3 <- dynsim.new(obj = int.mod.2, ldv = "l.vicab_filtered", scen = Scen3, n = 24,
                   shocks = shock_scen3)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)
Sim3<-as.data.frame(Sim3)


Sim1[5,4]-Sim1[11,4]
Sim2[5,4]-Sim2[11,4]
Sim3[5,4]-Sim3[11,4]

Sim1[5,4]
Sim1[11,4]
Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.113",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.16, yend=-0.95), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 1st quantile of", phantom(x), Delta, "unemployment")),x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]
Sim2[11,4]
Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.417",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.3, yend=-1.1), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the median of", phantom(x), Delta, "unemployment")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p


Sim3[5,4]
Sim3[11,4]
Sim3.p <- ggplot(Sim3, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-2,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.833",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.48, yend=-1.34), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 3rd quantile of", phantom(x), Delta, "unemployment")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim3.p

unemp.plot<-plot_grid(NULL,NULL,ins.dunemployment.austerity,NULL,NULL,NULL,
                      Sim1.p, NULL,
                      Sim2.p, NULL,
                      Sim3.p, NULL,
                      labels = c("","","A","","","",
                                 "B","","C","","D",""),label_size=20,nrow = 2, ncol = 6,align = 'h',rel_widths = c(2,0,2,0,2,0))
ggsave(filename="f1.eps", plot=unemp.plot,width = 28, height = 11,device=cairo_ps)



### Figure 2 austerity_6m * imf ###
attach(cabinet_filter_data)
set.seed(1234567)

int.mod.5<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_6m*imf,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.5)


## Instantaneous effect ##

pred_effects_by_regchange <- function(austerity_6m = 1){
  betas <- coef(int.mod.5)
  betas[1,2] <- 0.6291437
  betas<-as.numeric(betas[c(1), ])   
  varcov <- vcov(int.mod.5)
  beta.tilde <- data.frame(mvrnorm(1000,betas,varcov))
  beta.tilde <- beta.tilde[,c(2,3,8)]
  
  colnames(beta.tilde) <- c(
    "alpha_1","beta_0","beta_3") 
  
  
  
  # create storage for results
  effects <- data.frame(matrix(ncol=7,nrow=2))
  colnames(effects) <- c("inst.lo","inst.est","inst.hi","total.lo","total.est","total.hi","IMF_or_not")
  effects$IMF_or_not <- c("IMF","No_IMF")
  
  beta.tilde$sre_No_IMF  <- beta.tilde[,"beta_0"]
  
  beta.tilde$sre_IMF <- rowSums(beta.tilde[,c("beta_0",
                                              "beta_3")])
  
  effects$inst.lo[which(effects$IMF_or_not == "IMF")] <- unname(quantile(beta.tilde$sre_IMF, .025))
  effects$inst.est[which(effects$IMF_or_not == "IMF")] <- unname(summary(beta.tilde$sre_IMF)[4])
  effects$inst.hi[which(effects$IMF_or_not == "IMF")] <- unname(quantile(beta.tilde$sre_IMF, .975))
  effects$inst.lo[which(effects$IMF_or_not == "No_IMF")] <- unname(quantile(beta.tilde$sre_No_IMF, .025))
  effects$inst.est[which(effects$IMF_or_not == "No_IMF")] <- unname(summary(beta.tilde$sre_No_IMF)[4])
  effects$inst.hi[which(effects$IMF_or_not == "No_IMF")] <- unname(quantile(beta.tilde$sre_No_IMF, .975))
  
  return(effects)
}

effects <- pred_effects_by_regchange()
effects$regcode <- c(1,0)

require(dplyr)
distribution.dta<-select(cabinet_filter_data,c(d.unemployment,protest_freq,imf))



ins.imf.austerity <- ggplot(effects) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_segment(aes(x=regcode,xend=regcode,y=inst.hi,yend=inst.lo), size=1,show.legend=F) +
  geom_point(aes(x=regcode,y=inst.est,size=2.2),show.legend=F) +
  labs(title="Instantaneous effect of austerity_6m",x="",y="Vote Intention (cabinet)") + 
  geom_hline(yintercept=0, linetype=2) +
  scale_x_continuous(breaks=c(0,1), labels=c("No external involvement","External involvement"),limits=c(-.5,1.5)) + 
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20)) +
  geom_histogram (data=distribution.dta,aes(x = imf,y = (-..count..)*0.001),
                  inherit.aes = F, alpha = .3, fill="grey",
                  position = "stack", bins=30,binwidth = 0.1)
  

ins.imf.austerity

## Temporal illustration ## 
cabinet_filter_data[,"austerity_imf"]<-cabinet_filter_data$austerity_6m*cabinet_filter_data$imf
attach(cabinet_filter_data)
set.seed(1234567)
int.mod.5<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_imf,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.5)

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(cabinet_filter_data$protest_freq, na.rm = T),
                    austerity_imf=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=mean(cabinet_filter_data$protest_freq, na.rm = T),
                    austerity_imf=0)


shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_imf<- 0
shock_scen1$austerity_imf[6:11]<-0

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_imf<- 0
shock_scen2$austerity_imf[6:11]<-1


set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.5, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.5, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)



Sim1[5,4]-Sim1[11,4]
Sim2[5,4]-Sim2[11,4]

Sim1[5,4]
Sim1[11,4]

Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-4,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.357",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.25, yend=-1.09), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response without external involvement")),
       x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]
Sim2[11,4]
Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-4,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -1.8, label = "drop=3.680",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.259, yend=-3.42), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response with external involvement")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p



imf.plot<-plot_grid(ins.imf.austerity,NULL,NULL,NULL,
                    Sim1.p, NULL,
                    Sim2.p, NULL,
                    labels = c("A","","","","B","","C",""),label_size=20,nrow = 2, ncol = 4,align = 'h',rel_widths = c(2,0,2,0))

ggsave(filename="f2.eps", plot=imf.plot,width = 20, height = 12,device=cairo_ps)



### Figure 3 austerity_6m * protest ###
attach(cabinet_filter_data)
set.seed(1234567)

int.mod.8<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_6m*protest_freq,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.8)


## Instantaneous effect ##

pred_effects_d.protest_freq <- function(d.protest_freq){
  betas <- coef(int.mod.8)
  betas[1,2] <- 0.6284707
  betas<-as.numeric(betas[c(1), ])   
  varcov <- vcov(int.mod.8)
  beta.tilde <- data.frame(mvrnorm(1000,betas,varcov))
  beta.tilde <- beta.tilde[,c(2,3,8)]
  
  colnames(beta.tilde) <- c(
    "alpha_1","beta_0","beta_1") 
  
  # create relevant comparisons 
  protest_freq <- seq(unname(summary(cabinet_filter_data$protest_freq)[1]),
                      unname(summary(cabinet_filter_data$protest_freq)[6]),length.out=250) # vary level
  # create hypothetical values based on these
  hyp_zt <- protest_freq
  hyp_zt_l1 <- protest_freq - d.protest_freq
  hyp_zt_p1 <- protest_freq + d.protest_freq
  # create list to loop over
  hyplist <- as.list(NULL)
  # create storage for results
  effects <- data.frame(matrix(ncol=7,nrow=length(protest_freq)))
  colnames(effects) <- c("protest_freq","inst.lo","inst.est","inst.hi","total.lo","total.est","total.hi")
  for(i in 1:length(protest_freq)){
    ### create temp storage
    hyplist[[i]] <- beta.tilde
    
    ### create hypotheticals
    # lagged level
    
    hyplist[[i]]$beta_1_hyp_zt_l1 <- hyplist[[i]]$beta_1*hyp_zt_l1[i]
    
    # current level
    hyplist[[i]]$beta_1_hyp_zt <- hyplist[[i]]$beta_1*hyp_zt[i]

    # future level
    hyplist[[i]]$beta_1_hyp_zt_p1 <- hyplist[[i]]$beta_1*hyp_zt_p1[i]
    
    # negative future level
    hyplist[[i]]$beta_1_hyp_zt_p1_neg <- (-1)*hyplist[[i]]$beta_1*hyp_zt_p1[i]

    ### calculate quantities of interest
    hyplist[[i]]$inst.effect <- rowSums(hyplist[[i]][,c("beta_0", "beta_1_hyp_zt" )])
    
    hyplist[[i]]$total.effect <- (rowSums(hyplist[[i]][,c("beta_0", "beta_1_hyp_zt"
    )]))/(1-(hyplist[[i]]$alpha_1))
    effects$protest_freq[i] <- hyp_zt[i]
    effects$inst.lo[i] <- unname(quantile(hyplist[[i]]$inst.effect, .025))
    effects$inst.est[i] <- unname(summary(hyplist[[i]]$inst.effect)[4])
    effects$inst.hi[i] <- unname(quantile(hyplist[[i]]$inst.effect, .975))
    effects$total.lo[i] <- unname(quantile(hyplist[[i]]$total.effect, .025))
    effects$total.est[i] <- unname(summary(hyplist[[i]]$total.effect)[4])
    effects$total.hi[i] <- unname(quantile(hyplist[[i]]$total.effect, .975))
  }
  return(effects)
}

effects <- pred_effects_d.protest_freq(d.protest_freq = 0)

require(dplyr)
distribution.dta<-select(cabinet_filter_data,c(d.unemployment,protest_freq,imf))

ins.protest.austerity <- ggplot(effects, aes(x=protest_freq,y=inst.est)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=inst.lo, ymax=inst.hi), alpha=.15) +
  geom_line(aes(x=protest_freq,y=inst.est)) +
  coord_cartesian(ylim=c(5,2.5), xlim = c(0,3)) + 
  geom_hline(yintercept=0, linetype=2) + 
  labs(title="Instantaneous effect of austerity_6m", x=expression("Protest Frequency"),y="Vote Intention (cabinet)") + 
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20)) +
  geom_rug (data=distribution.dta,aes(x = protest_freq),
            inherit.aes = F, alpha = .5,
            position = "identity",size=1,sides="b")
ins.protest.austerity



ins.protest.austerity <- ggplot(effects, aes(x=protest_freq,y=inst.est)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=inst.lo, ymax=inst.hi), alpha=.15) +
  geom_line(aes(x=protest_freq,y=inst.est)) +
  coord_cartesian(ylim=c(-5,2.5), xlim = c(0,3)) + 
  geom_hline(yintercept=0, linetype=2) + 
  labs(title="Instantaneous effect of austerity_6m", x=expression("Protest Frequency"),y="Vote Intention (cabinet)") + 
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20)) +
  geom_histogram (data=distribution.dta,aes(x = protest_freq,y = (-..count..)*0.005),
                  inherit.aes = F, alpha = .3, fill="grey",
                  position = "stack", bins=30,binwidth = 0.1)
ins.protest.austerity


## Temporal illustration ## 

cabinet_filter_data[,"austerity_protest_freq"]<-cabinet_filter_data$austerity_6m*cabinet_filter_data$protest_freq
attach(cabinet_filter_data)
set.seed(1234567)
int.mod.8<- lme(
  fixed = vicab_filtered ~ l.vicab_filtered+austerity_6m
  + d.unemployment+ retail+ imf+protest_freq+austerity_protest_freq,
  random = ~ 1+l.vicab_filtered | country, control = lmeControl(opt = 'optim'))
summary(int.mod.8)

Scen1 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.1, na.rm = T),
                    austerity_protest_freq=0)

Scen2 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.75, na.rm = T),
                    austerity_protest_freq=0)

Scen3 <- data.frame(l.vicab_filtered = 0,
                    austerity_6m=0,
                    d.unemployment = mean(cabinet_filter_data$d.unemployment, na.rm = T),
                    retail = mean(retail,
                                  na.rm = TRUE),
                    imf=0,
                    protest_freq=quantile(cabinet_filter_data$protest_freq,probs =0.9, na.rm = T),
                    austerity_protest_freq=0)

shock_scen1 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen1$austerity_6m [6:11]<- 1 
shock_scen1$austerity_protest_freq<- 0
shock_scen1$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.1, na.rm = T)

shock_scen2 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen2$austerity_6m [6:11]<- 1 
shock_scen2$austerity_protest_freq<- 0
shock_scen2$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.75, na.rm = T)


shock_scen3 <- data.frame(times = 1:24, austerity_6m = 0)
shock_scen3$austerity_6m [6:11]<- 1 
shock_scen3$austerity_protest_freq<- 0
shock_scen3$austerity_protest_freq[6:11]<-1*quantile(cabinet_filter_data$protest_freq,probs =0.9, na.rm = T)

set.seed(1234567)
Sim1 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen1, n = 24,
                   shocks = shock_scen1)
set.seed(1234567)
Sim2 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen2, n = 24,
                   shocks = shock_scen2)
set.seed(1234567)
Sim3 <- dynsim.new(obj = int.mod.8, ldv = "l.vicab_filtered", scen = Scen3, n = 24,
                   shocks = shock_scen3)
Sim1<-as.data.frame(Sim1)
Sim2<-as.data.frame(Sim2)
Sim3<-as.data.frame(Sim3)


Sim1[5,4]-Sim1[11,4]
Sim2[5,4]-Sim2[11,4]
Sim3[5,4]-Sim3[11,4]

Sim1[5,4]
Sim1[11,4]
Sim1[5,4]-Sim1[11,4]
Sim1.p <- ggplot(Sim1, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.518",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.49, yend=-1), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 10% quantile of protest frequency")),x="Month",y=expression(paste("Predicted Vote Intention (Cabinet)"))) +
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim1.p

Sim2[5,4]
Sim2[11,4]
Sim2[5,4]-Sim2[11,4]
Sim2.p <- ggplot(Sim2, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=1.648",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.41, yend=-1.23), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 75% quantile of protest frequency")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim2.p


Sim3[5,4]
Sim3[11,4]
Sim3[5,4]-Sim3[11,4]
Sim3.p <- ggplot(Sim3, aes(x=time,y=ldvMean)) + 
  theme_bw(base_size =24) + 
  theme(axis.line = element_line(colour = "black"),panel.border = element_blank(),panel.grid.minor = element_blank()) +
  geom_ribbon(aes(ymin=ldvLower, ymax=ldvUpper), alpha=.15) +
  geom_ribbon(aes(ymin = ldvLower50, 
                  ymax = ldvUpper50), alpha = .1, linetype = 0) +
  geom_line(aes(x=time,y=ldvMean)) + 
  scale_x_continuous(limits=c(1,24),
                     breaks=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24),
                     labels=c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24)) +
  coord_cartesian(ylim=c(-3,1),xlim=c(1,24)) + 
  geom_hline(yintercept=0, linetype=2) +
  annotate("text", x = 3, y = -0.6, label = "drop=2.145",size = 6) +
  geom_segment(inherit.aes = F, aes(x=5, xend=5, y=0.11, yend=-2.03), 
               size=0.3, linetype=2, arrow = arrow(length = unit(0.1, "cm"),ends = "both", type = "closed")) +
  labs(title=expression(paste("Impulse response at the 90% quantile of protest frequency")),x="Month",y="")+
  theme(plot.title = element_text(size=24),axis.title= element_text(size=20))

Sim3.p


protest.plot<-plot_grid(NULL,NULL,ins.protest.austerity,NULL,NULL,NULL,
                        Sim1.p, NULL,
                        Sim2.p, NULL,
                        Sim3.p, NULL,
                        labels = c("","","A","","","",
                                   "B","","C","","D",""),label_size=20,nrow = 2, ncol = 6,align = 'h',rel_widths = c(2,0,2,0,2,0))

ggsave(filename="f3.eps", plot=protest.plot,width = 30, height = 12,device=cairo_ps)




